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We introduce a new linear response method to study the lattice dynamics of materials with strong 
correlations. It is based on a combination of dynamical mean field theory of strongly correlated 
electrons and the local density functional theory of electronic structure of solids. We apply the 
method to study the phonon dispersions of a prototype Mott insulator NiO. Our results show 
overall much better agreement with experiment than the corresponding local density predictions. 
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Computational studies of lattice dynamics and struc- 
tural stability in strongly correlated electronic systems 
is a challenging theoretical problem. In the past den- 
sity functional theory (DFT) in its local density or 
generalized gradient approximations (LDA or GGA) Jl]] 
has delivered the full lattice dynamical information and 
electron-phonon related properties of a variety of sim- 
ple metals, transition metals, as well as semiconductors 
with exceptional accuracy 0] . This is mainly due to the 
introduction of a linear response approach This 
method overcame the problems of traditional techniques 
based on static susceptibility calculations which gener- 
ally fail to reproduce lattice dynamical properties of real 
materials due to difficulties connected with the summa- 
tions in high-energy states and the inversion of very large 
dielectric matrix ||. 

Despite these impressive successes, there is by now 
clear evidence that the present methodology fails when 
applied to strongly correlated materials. For example, 
the local density predictions for such properties as bulk 
modulus and elastic constants in metallic Plutonium are 
approximately one order of magnitude off from experi- 
ment H; the phonon spectrum of Mott Insulators such 
as MnO is not predicted correctly by LDA . 

In this work we describe a new linear response method 
to study the lattice dynamics of correlated materials. It 
is based on the dynamical mean field theory (DMFT) ||, 
a many-body technique developed to study systems with 
strong on-site Coulomb repulsion. Recent progress in 
merging |J this many-body description with the realistic 
LDA based electronic structure calculations has already 
led to new insights in long standing problems, such as, 
the temperature dependence of the magnetic properties 
of Fe and Ni [jHJ, the volume collapse transition in Ce 
0, and Pu ifQ. We generalize this LDA+DMFT 
method to carry out linear response calculations by find- 
ing self-consistent changes in both charge densities and 
local Green functions induced by atomic displacements. 

The LDA+DMFT approach computes the total energy, 
the charge density and the local spectral function of the 



correlated electrons simultaneously. The latter spectra 
are known to be quite different for a strongly correlated 
system such, e.g., as heavy fermion metal, from the Kohn 
Sham spectra of the local density functional theory due 
to the appearance of strongly renormalized quasiparticle 
features at the Fermi level and lower and upper Hub- 
bard features or satellites at higher energies ||. The 
LDA+DMFT technique can therefore provide a link be- 
tween the photoemission spectra and the lattice dynam- 
ics of correlated materials, which is an interesting open 
problem. 

As a test we consider the prototype Mott insulator 
NiO. In this material at low temperatures, the static 
Hartree-Fock limit of our formulation which is formally 
equivalent to the LDA+U method E3] can be used. Re- 
sults for the phonon dispersion and the static dielec- 
tric properties of this material are found in much better 
agreement with the experiment than the results of corre- 
sponding LDA based calculations. We discuss how corre- 
lations affect the calculated dielectric constants and the 
Born effective charges using the linear-response method. 

Our approach considers both the charge density p and 
the local Green function G(u) as parameters of a spectral 
density functional fL3| . To find its extremum, a set of 
Dyson equations is solved self-consistently: 



(-V 2 + V ( 



eff 







(1) 



where V e ff is the effective potential of the local den- 
sity functional and E(w) is the local self-energy oper- 
ator. Since the local density approximation contains an 
average correlation energy, a double counting term T,dc 
appears in the equation (|l|). Both V e // and are 
functionals of the density and the local Green function 
which can be found using the formula: 



(2) 



Note, that the Greens function is a non-hermitian ma- 
trix that is frequency dependent so both the eigenvalues 
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e and right/left eigenvectors tp r ' are treated formally as 
frequency dependent quantities: ^kju)i' t Pkj u) , ^kju- (The 
latter satisfies the Dyson equation ([!]) with the wave func- 
tion placed on the left.) In practice ||Jl^], (0) ^ s solved 
on the Matsubara axis for a finite set of imaginary fre- 
quencies ioj n using some localized orbital representation 
such, e.g., as linear muffin-tin orbitals (LMTOs) (lBJXa 
for the eigenvectors ip^^ 



(5) 
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which substitutes the differential equation ([!]) by a ma- 
trix eigenvalue problem. 

Once the local Green function is constructed, the new 
charge density, effective potential and the local self- 
energy are computed. The latter is found by solving 
the Anderson impurity model using a suitable many- 
body technique. This entire formulation requires the self- 
consistent procedure which delivers the total energy of 
the interacting electronic system. 

The dynamical matrix is the second order derivative 
of the energy. As with the ordinary density functional 
formulation of the problem Jl6|], we deal with the first 
order corrections to the charge density, 8p, as well as the 
first order correction to the local Greens function SG(ui) 
which should be considered as two independent variables 
in the functional of the dynamical matrix. To find the 
extremum, a set of the linearized Dyson equations has to 
be solved self-consistently: 



(-V 2 + V e ff + EM - E dc - e kj Ml ju + 

+(sv eff + s£(u) - s£ dc )ipi jul = o 



(4) 



which leads us to consider the first order changes in the 
effective potential 5 V e f / and in the local self-energy op- 
erator S~E(lo). Here and in the following we will assume 
that the phonon wave vector of the perturbation q is dif- 
ferent from zero, and, therefore, the first order changes 
in the eigenvalues Sc^joj drop out. The quantities SV e ff 
and <$E(w) are the functionals of 5p and 6G(lu) and should 
be found self-consistently. In particular, the change in 
the self-energy #E(a)) assumes the development of solv- 
ing an Anderson impurity model linearized with respect 
to atomic displacements. 

In practice, change in the eigenvector 8ipkj U has to 
be expanded in some basis set. Previous linear re- 
sponse schemes were based on tight-binding methods 
p7| , plane wave pseudopotentials P^, p^ , p^| , linear aug- 
mented plane wave s [pp| , mixed orbitals ]2l| and linear 
muffin-tin orbitals ]22|. To build an effective computa- 
tional scheme applicable for systems with localized or- 
bitals we use LMTO representation as the basis. Due 
to its explicit dependence on the atomic positions both 
Hellmann-Feynman contributions and incomplete basis 
set corrections [ p3| appear in the expression for the dy- 
namical matrix ]16fl. We expand Sip^j^ as follows 



where we introduced both changes in the frequency de- 
pendent variational coefficients 5A^ U as well as changes 
in the basis functions #Xa- The latter helps us to reach 
convergency in the entire expression (^) with respect to 
the number of the basis functions {a} fast since the con- 
tribution with (5%^ takes into account all rigid movements 
of the localized orbitals . 

The first-order changes in the Green function can be 
found as follows 



(6) 



which should be used to evaluate the first order change 
in the charge density and the dynamical matrix itself. 

Despite formal complexity of the outlined method we 
see two significant simplifications which can be imple- 
mented in practice. First, using static limit of the dynam- 
ical mean field theory leads to a frequency independent 
approximation for E. This recovers the LDA+U method 
and its linear-response analog for the present problem. 
Second, it is generally expected that only low frequency 
behavior of the self-energy will drastically affect the static 
electronic response, since it modifies the one-electron 
spectrum near the Fermi level. If the Hubbard subbands 
are separated by sufficiently large U they are essentially 
quasi-atomic features and may be though as the states 
tightly bound and moving rigidly as the atom vibrates. 
These two simplified calculational schemes may be ex- 
plored for the lattice dynamics in correlated materials 
without significant complications. 

We now describe our implementation of the method for 
calculating vibrational spectrum of the prototype Mott 
insulator NiO. It is well known |24j] that the LDA un- 
derestimates both the value of the energy gap and the 
magnetic moment of this material. On the other hand, 
the use of the LDA+U method fixes both problems |2~i| ] 
when using non-zero Hubbard parameter U. The use of 
the LDA+U approximation is equivalent here to a static 
limit of the full frequency resolved theory outline above: 
the self-energy E(w) becomes an orbital-dependent cor- 
rection to the Kohn-Sham potential U e // expressed via 
the density matrix of the localized electrons, n Q/ g, and the 
change in the self-energy, <5E, is expressed via the changes 
in occupation numbers, Sn a /3, due to atomic movements. 
Here we will assume that the on-site Hubbard interac- 
tion U is fixed for all the displacements. The use of the 
static approximation significantly simplifies our imple- 
mentation: the solution of the impurity model and its 
change with respect to the displacements collapses and 
the quantities which have to be found self-consistently are 
the changes in charge density bp and in the occupancy 
matrix 5n a j3. The latter substitutes finding full change 
SG(u>) in the local Green function of the system. 
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FIG. 2. Comparison between calculated using the present 
method (filled circles) and experimental (open circles) phonon 
dispersion curves for NiO. 

We first discuss the results of the LDA calculations 
for the lattice dynamics of NiO. Fig. 1 shows calculated 
phonon dispersion relations along major symmetry di- 
rections. Despite the unit cell corresponding to antifer- 
romagnetic lattice contains two original unit cells, spin 
integrated properties are unaffected by this symmetry 
lowering and therefore can be represented in the original 
Brillouin zone. We see an apparent underestimation of 
all theoretical optical modes by as much as 50 or more 
per cent as compared to the experiment pq|. Similar 



finding for MgO has also been reported fjj . The LO-TO 
splitting is largely underestimated which is predicted to 
be only 1.5 THz against 5 THz seen from the measure- 
ment. Surprisingly, however, that both longitudinal and 
transverse acoustic branches are much better predicted 
by the LDA. One can conclude that the calculated sound 
velocities along TX direction are only slightly overesti- 
mated and the agreement is even better along the TL 
line. 

We now discuss how inclusion of correlations modifies 
the phonon spectrum. Fig. 2 shows the calculated dis- 
persion curves using static self-energy approximation of 
the formalism developed in this work. We use the val- 
ues of Hubbard U=8 eV and exchange J=l eV in this 
calculation. We see a remarkable hardening of the opti- 
cal modes which occurs due to proper treatment of the 
correlational effects. Notably both the values of the op- 
tical frequencies as well as the calculated LO-TO split- 
ting are now in a good agreement with the experiment. 
A pronounced softening of the longitudinal optical mode 
along both TX and TL lines is seen at the measured data 
which is in part captured by our theoretical calculation: 
the agreement is somewhat better along the TX direction 
while the detailed q-dependence of these branches shows 
some residual discrepancies. 

Our acoustic modes are consistently overestimated by 
a few per cent as compared to the experiment. They 
are also harden as the on-site interaction U increases al- 
though in much smaller degree than the optical modes. 
Overall, the agreement between the theory and the ex- 
periment is much better when the correlations are taken 
into account which indicates their importance in the lat- 
tice dynamics of strongly correlated materials. 

To better understand our findings we have calculated 
the values of which represent an electronic contri- 
bution to the static dielectric constant as given by the 
inverse element l/e _1 (0,0) of the full inverse dielectric 
permittivity matrix in reciprocal space. Table 1 shows 
the results of these calculations using the LDA and our 
new method as well as the experiment pffl . It is seen that 
the LDA result overestimates the value of by a fac- 
tor of 6. This is mainly because the calculated by LDA 
dielectric gap value is only 0.4 eV which is much smaller 
than the experimental energy gap of NiO equal to 4 eV 
p8| . Namely LDA overestimates the electronic screening 
effects by a large amount causing both the artificial soft- 
ening of optical phonons and the lowering of the LO-TO 
splitting. The latter is well known to be directly pro- 
portional to the Born effective charges but inversionally 
proportional to e,^. 

On the other hand, our calculations with correlations 
produce a much better value of the static dielectric con- 
stant as seen from Table 1. We relate such an agreement 
with the fact that the energy gap value predicted by the 
theory here is 3.3 eV which compares well with the mea- 



surement 28 1 . This indicates that the screening effects 
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in the calculation with correlations are treated in much 
more appropriate way which fixes both the position of 
the optical branches and their relative splitting. 

As a final result, we have extracted the values of the 
Born effective charge Z* using q — ► limit technique 



29| and the relationship lo 
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|Z*|7 eoo . Table 1 



shows comparison between the theoretical and the exper- 
imental pt| data. The deviation from 2 which for binary 
oxides is the nominal rigid-ion value of \Z* | indicates high 
electronic polarizability. Both theoretical values are close 
to the experimental one, while the agreement for LDA is 
apparently accidental since the error in underestimating 
the splitting ui\ — oj^o nearly completely cancels the 
error in overestimating eoo. 

To summarize we have developed a new method to 
study the lattice dynamical properties of materials with 
strong electronic correlations. The method is based on 
the dynamical mean-field theory and it uses linear re- 
sponse technique to calculate phonons at arbitrary wave 
vector q. Using static approximation we have computed 
the phonon dispersions of NiO and find the results for 
the optical modes and static dielectric properties in much 
better agreement with the experiment than the results of 
the corresponding LDA based calculations. This empha- 
sizes the importance of correlations in lattice dynamics 
of strongly correlated systems. The residual discrepan- 
cies can be attributed to the frequency dependence of 
the self-energy neglected in the present study and will 
be a subject for the future work. There are many other 
challenging problems for which the general framework 
outlined here can be useful. Discrepancies between the- 
ory and experiment were noticed p0[ | for the ferroelec- 
tric CaCu3Ti40i2- The electron-phonon interaction in 
cuprate superconductors is a subject of intensive inves- 
tigation. Finally we believe that correlation effects play 
an important role in the phonon dynamics across the ac- 
tinide series. 
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TABLE I. Comparison between calculated static dielectric 
constant and Born effective charge \Z*\ for NiO using LDA 
and present method as well as the experimental data |2q,p7f 





LDA 


Present Method 


Exp 


Coo 


35.7 


7.2 


5.7, 6.1 


\z*\ 


2.17 


2.33 


2.22 
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